library(terra)
library(lidR)
library(sf)
library(here)
library(data.table)
zmax <- 40
grain.size <- 1
project_root <- here::here()
# Choose LAD method: "linear" or "beer"
# Beer–Lambert conversion Notes:
# - Avoids log(0) and 1 by clipping near-extreme values
# - Use when cumulative light absorption or occlusion is relevant
# - Suitable if extinction coefficient is known or estimated from prior studies
lad_method <- "beer" # Set to "linear" or "beer"
# Optional: extinction coefficient (used only for Beer–Lambert conversion)
k_extinction <- 0.25
las_file <- file.path(project_root, "data/TLS/tree_08.laz")
output_voxels <- file.path(project_root, "data/TLS/LAD_voxDF.rds")
output_array <- file.path(project_root, "data/TLS/lad_array_m2m3.rds")
output_profile_plot <- file.path(project_root, "data/TLS/lad_vertical_profile.pdf") Using leaf area density (LAD) from TLS data in ENVI-met for 3D plants.
Background and Method
This section explains the theoretical principles of leaf area density (LAD) and describes how it can be determined using terrestrial laser scanning (TLS). Leaf area density is an important parameter in environmental modeling, for example for radiation balance and microclimate simulations. It indicates the leaf area per volume (m²/m³) and is therefore a decisive factor for microclimate simulations, radiation models, and energy flows in vegetation stands.
| Approach Type | Name / Description | Nature |
|---|---|---|
| Pulse-count based | Simple linear normalization of return counts or voxel hits | Empirical, direct |
| Linear normalization | Straightforward normalization of pulse counts by voxel volume or max LAD | Empirical, basic |
| Pulse-density normalization | Adjusts for occlusion and scan geometry | Semi-empirical |
| Gap fraction models | Estimate LAD/LAI from canopy openness statistics | Semi-empirical |
| Beer–Lambert conversion conversion | Uses exponential light attenuation to infer LAD | Physically-based |
| Voxel-based inverse modeling | Optimizes 3D LAD to match observed light attenuation or reflectance | Physically-based |
| Allometric / geometric reconstruction | Reconstructs crown volume and distributes LAD using QSM or shape fitting | Geometric, structural |
- Linear normalization is a practical baseline: simple, fast, and reproducible.
- Beer–Lambert conversion introduces realism via physical light attenuation.
More advanced models (e.g. voxel inverse or QSM-based) aim for higher biophysical fidelity at the cost of complexity.
The present analysis is based on TLS with a medium-range RIEGL scanner (e.g., VZ-400). This captures millions of 3D points of the vegetation structure with high angular resolution. The point cloud is divided into uniform voxels, from which the leaf area density is estimated in two ways.
Linear normalization (straightforwad)
\[
\text{LAD}_i = \frac{N_i}{N_{\max}} \cdot \text{LAD}_{\max}
\] - \(N_i\): Number of laser points in voxel \(i\)
- \(N_{\max}\): Maximum across all voxels
- \(\text{LAD}_{\max}\): Maximum LAD value from the literature (e.g., 5 m²/m³)
Beer–Lambert conversion
\[ \text{LAD}_i = -\frac{\ln\left(1 - \frac{N_i}{N_{\max}}\right)}{k \cdot \Delta z} \]
- \(k\): Extinction coefficient (typically 0.3–0.5)
- \(\Delta z\): vertical voxel height
Overall Workflow
What happens in the script?
| Step | Description | Relevant Code |
|---|---|---|
| 1. Read & Filter LAS | Load TLS data, optionally crop and clean it | readLAS() and las = filter_poi(...) |
| 2. Voxel Grid Setup | Set up 3D grid at defined grain.size |
passed to pixel_metrics(..., res = grain.size) |
| 3. Count Pulses | Count returns in each voxel height bin | pointsByZSlice() function |
| 4. Normalise Pulse Counts | Divide by global max (relative LAD) | in convert_to_LAD(): lad = (count / max) * LADmax |
| 5. Export Raster | Convert metrics to raster stack | terra::rast() from voxel_df |
| 6. Visualization | Plot LAD profiles | see plotting section |
| 7. Export to Plant3D | Exports the LAD to ENVI-met | see export section |
Implemetation
To use this ENVI-met tree modeling workflow in R, follow these steps to load and initialize the project correctly:
Project Setup: Loading the R Project and Environment
Download and Unzip the Project Archive
Unzip the folder to your desired location.
The folder should contain at least:
- An
*.Rprojfile (e.g.envimet_tree_workflow.Rproj) - A
data/folder with input files liketree_08.las - One or more
R/scripts
- An
Open the Project in RStudio
- Go to File → Open Project
- Select the
*.Rprojfile (e.g.microclimate_TLS.Rproj) - This ensures that the project directory is treated as the root for all file paths.
The use of the
{here}package depends on having a valid RStudio project. Without this, file paths may not resolve correctly.
Data Input Parameters and Paths
Set global parameters for the workflow, such as file paths, voxel resolution, and maximum LAD value for normalization.
Voxelization of TLS data
Voxelisation turns a 3D TLS point cloud into a grid of cubes (voxels), where each voxel holds structural information. The number of points per voxel is used to estimate Leaf Area Density (LAD), typically normalized relative to the voxel with the most returns.
- Each voxel = a 1×1×1 m³ cube
- Count the laser hits per voxel
- Normalize to maximum
- Multiply by a literature-based LAD_max (e.g. 5 m²/m³)
This gives a spatially distributed LAD profile suitable for further analysis or models like ENVI-met.
Conversion to LAD (m²/m³)
The conversion to LAD (Leaf Area Density, in m²/m³) from TLS-based voxel pulse counts is done using a relative normalization heuristic which is adopted as a practical approximation in voxel-based canopy structure analysis using TLS (Terrestrial Laser Scanning) data.:
For each voxel layer (e.g. pulses_2_3m), the LAD is calculated as:
\[ \text{LAD}_{\text{voxel}} = \left( \frac{\text{pulse count in voxel}}{\text{maximum pulse count over all voxels}} \right) \times \text{LAD}_{\text{max}} \]
Where:
pulse count in voxel= number of returns in this voxel layer (from TLS)max_pulse= the maximum pulse count found in any voxel (used for normalization)LAD_max= a fixed normalization constant (e.g. 5.0 m²/m³) chosen from literature or calibration
Typical LADₘₐₓ Values by Species
| Species / Structure Type | LADₘₐₓ (m²/m³) | Source / Notes |
|---|---|---|
| Fagus sylvatica (European beech) | 3.5–5.5 | Calders et al. (2015), Chen et al. (2018) |
| Quercus robur (English oak) | 3.0–6.0 | Hosoi & Omasa (2006), field studies with TLS voxelization |
| Coniferous trees (e.g. pine) | 4.0–7.0 | Wilkes et al. (2017), higher LAD due to needle density |
| Mixed broadleaf forest | 3.0–6.0 | Flynn et al. (2023), canopy averaged estimates |
| Shrubs / understorey | 1.5–3.0 | Chen et al. (2018),lower vertical structure density |
| Urban street trees | 2.0–4.0 | Simon et al. (2020), depending on pruning and species |
LAD values refer to maximum expected per 1 m vertical voxel. Values depend on species, seasonality, and scanning conditions.
What this means conceptually
You’re not measuring absolute LAD, but instead:
- Using the number of TLS returns per voxel as a proxy for leaf density
- Then normalization all voxels relatively to the most “leaf-dense” voxel
- The
LAD_maxdefines what value the “densest” voxel should reach in terms of LAD
This is fast, simple, and works well when:
- You want relative structure across the canopy
- You don’t have absolute calibration (e.g. with destructive sampling or hemispheric photos)
Caveats and assumptions
- This approach assumes the TLS beam returns are proportional to leaf area, which is a simplification
- It’s sensitive to occlusion and TLS positioning
- The choice of
LAD_maxis crucial—common values from literature range from 3–7 m²/m³ for dense canopies
The LAD conversion in the following code is a relative, normalized mapping of TLS pulse counts to LAD values, normalized by the highest voxel return and normalized using a fixed LAD_max. This gives a plausible LAD field usable for analysis, visualization, or simulation input (e.g. for ENVI-met).
DT::datatable(head(df_lad, 5))Raster Stack Representation of 3D Vegetation (Voxel-Based)
We represent 3D vegetation using a voxel-based raster stack:
- Space is divided into cubic voxels (e.g. 1 × 1 × 1 m).
- Each raster layer represents a height slice (e.g. 0–1 m, 1–2 m, …).
- Voxels store values like pulse counts or Leaf Area Density (LAD).
This 2D stack structure enables:
- Vertical profiling of vegetation per XY column.
- Layer-wise analysis (e.g. median, entropy).
- Integration with raster data like topography or irradiance.
- Use in raster-based ecological and microclimate models.
It supports both analysis and visualization of vertical structure with standard geospatial tools.
ENVI-met supports custom vegetation input via the SimplePlant method, which requires a vertical LAD profile per grid column. A raster stack derived from TLS data provides exactly this: each layer represents LAD in a specific height slice, and each XY cell corresponds to one vertical profile. This structure can be exported as CSV, ASCII rasters, or custom profile files.
For 3D vegetation parameterization in ENVI-met 5.8+, the raster stack enables preprocessing of spatially explicit LAD or LAI profiles, even if some reformatting is needed.
The raster stack also supports canopy clustering and prototyping. It allows classification of structural types, simplification of complex vegetation, and the creation of representative profiles for simulation.
# In SpatRasterStack umwandeln
xy <- df_lad[, c("X", "Y")]
lad_vals <- df_lad[, grep("^lad_", names(df_lad), value = TRUE)]
lad_raster <- rast(cbind(xy, lad_vals), type = "xyz")
plot(lad_raster)
In a more 3D version it looks like below.
null
1 Visualization
LAD Profile Visualizations from TLS Data
The plot_lad_profiles() function visualizes vertical leaf area density (LAD) profiles derived from voxelized TLS (terrestrial laser scanning) data. LAD represents leaf surface area per unit volume (m²/m³). The function provides three main plot styles:
1. XY Matrix Plot (plotstyle = "each_median")
Displays a grid of mini-profiles, each representing a 0.5 × 0.5 m (x/y) ground column.
Within each cell, a normalized vertical LAD profile is plotted:
- Y-axis (height) is normalized from 0 to 1 per column.
- X-axis shows LAD values normalized relative to the global LAD maximum.
Useful for comparing structural patterns across space.
2. Overall Median Profile (plotstyle = "all_median")
- Aggregates LAD values across all (x/y) locations by height bin.
- Produces a typical vertical profile using the median and smoothed with a moving average.
- Height is shown in absolute units (e.g. meters).
- Captures the dominant vertical canopy structure.
3. Single Profile (plotstyle = "single_profile")
- Extracts and plots the LAD profile at a specific (x, y) coordinate.
- Both LAD and height are shown in absolute units.
- Plots the true vertical structure at one location.
The matrix plot shows multiple vertical LAD profiles arranged in a grid, with each small plot corresponding to a specific spatial location. This allows the vertical vegetation structure to be viewed in relation to its position on the ground. To make the individual profiles comparable, both height and LAD values are normalized within the plot. A reference profile on the side shows the overall median LAD distribution by height, which helps interpret the scale and shape of the individual profiles.
# Option 1: Profile in each column
plot_lad_profiles(lad_df, plotstyle = "each_median")
# Option 2: Overall vertical LAD profile (median of all)
plot_lad_profiles(lad_df, plotstyle = "all_median")
# Option 3: Single profile at specified coordinates
plot_lad_profiles(lad_df, plotstyle = "single_profile", single_coords = c(57.5, -94.5))
ENVI-met 3D Tree Export
The next section describes more detailed how the key input values in the R function export_lad_to_envimet3d() are computed, derived or selected, and provides the rationale for each. The function converts a voxel-based Leaf Area Density (LAD) profile, typically obtained from Terrestrial Laser Scanning (TLS) data, into a structured XML file compatible with ENVI-met’s 3D tree model (.pld or PLANT3D).
Given the sensitivity of ENVI-met simulations to tree morphology and LAD distribution, the function ensures that the spatial dimensions, vertical layering and LAD intensity values are all correctly represented. Some parameters are optional, but can be derived from the data if not explicitly set.
The table below details each argument of the function, including its purpose, how it is determined and its necessity.
| Code Line | Meaning | Reason |
|---|---|---|
lad_df <-lad_df[!is.na(lad_df$LAD_median), ] |
Removes entries with missing LAD values | Ensures only valid data is used in the LAD calculation and XML export |
lad_df$i <-as.integer(factor(lad_df$x)) |
Converts x-coordinates to integer voxel column indices (i) | Required for ENVI-met LAD matrix indexing |
lad_df$j <-as.integer(factor(lad_df$y)) |
Converts y-coordinates to integer voxel row indices (j) | Same as above, for the y-direction |
z_map <-setNames( ...) |
Maps unique height bins to sequential vertical indices (k) | Translates height levels into voxel layers compatible with ENVI-met |
lad_df$k <-z_map[as.character(lad_df$Height_bin)] |
Applies the vertical index to the LAD data | Aligns LAD values with ENVI-met vertical layer system |
lad_df$lad_value <-round(lad_df$LAD_median * scale_factor, 5) |
Scales LAD values and rounds to 5 digits | Brings LAD values to a usable range for ENVI-met and ensures precision |
dataI <-max(lad_df$i) |
Gets the number of horizontal grid cells in i-direction (width) | Required as matrix size input for ENVI-met |
dataJ <-max(lad_df$j) |
Gets the number of horizontal grid cells in j-direction (depth) | Required as matrix size input for ENVI-met |
zlayers <-max(lad_df$k) |
Gets the number of vertical layers | Sets the height resolution of the LAD matrix |
Automatic Grid Dimensions transformation
Calculates the voxel grid dimensions in X, Y, and Z from the TLS-derived LAD profile.
The table below outlines how the core spatial and structural parameters of the tree model are computed from the input LAD_DF data frame. These derived values define the three-dimensional structure of the tree in terms of its horizontal extent, vertical layering and canopy dimensions.
Data I and data J represent the size of the voxel grid in the i and j dimensions, respectively, based on unique horizontal (x and y) and vertical (height bin) bins in the LAD profile.
‘Width’ and ‘Depth’ describe the physical spread of the tree crown, inferred from the voxel grid extent if not manually set.
Height is computed by multiplying the number of vertical layers (zlayers) by the voxel resolution (cellSize), providing the total modelled height of the canopy.
These computed values are essential for correctly normalization and locating the 3D LAD matrix within the ENVI-met simulation domain to ensure visual and physiological realism.
| Code Line | Meaning | Reason |
|---|---|---|
Width <- if (is.null(Width)) dataI else Width |
Uses the number of i-cells if Width is not provided | Automatically estimates tree width from voxel spread in x-direction |
Depth <- if (is.null(Depth)) dataJ else Depth |
Uses the number of j-cells if Depth is not provided | Automatically estimates tree depth from voxel spread in y-direction |
Height <- zlayers * cellsize |
Converts number of vertical layers to metric height using cellsize | Computes physical tree height in meters for ENVI-met |
# 1. Remove NA values from the LAD column
lad_df <- lad_df[!is.na(lad_df$LAD_median), ]
# 2. Create discrete i and j indices for the horizontal position
# (converts x and y coordinates into consecutive index values)
lad_df$i <- as.integer(factor(lad_df$x))
lad_df$j <- as.integer(factor(lad_df$y))
# 3. Assign each Height_bin (z direction) a consecutive layer ID k
# (z_map assigns an index layer to each unique height)
z_map <- setNames(seq_along(sort(unique(lad_df$Height_bin))), sort(unique(lad_df$Height_bin)))
lad_df$k <- z_map[as.character(lad_df$Height_bin)]
# 4. Scale LAD values, e.g. to get from 0.02 to more realistic values such as 0.5–1.5
lad_df$lad_value <- round(lad_df$LAD_median * scale_factor, 5)
# 5. Calculate the maximum dimensions of the grid (for XML specifications)
dataI <- max(lad_df$i) # Width in cells (x-direction)
dataJ <- max(lad_df$j) # Depth in cells (y-direction)
zlayers <- max(lad_df$k) # Number of vertical layers (z-direction)Transmittance and Albedo
Albedo = 0.18
Transmittance = 0.3Albedo = 0.18: Albedo is the fraction of incoming solar radiation reflected by the canopy surface. For deciduous trees, values usually range between 0.15 and 0.20. 0.18 is a commonly used default for broadleaved species like Fagus sylvatica or Quercus robur in many ecological models (e.g., ENVI-met, MAESPA). It affects surface energy balance and radiation reflection in ENVI-met simulations.
Transmittance = 0.3: Transmittance represents the proportion of shortwave radiation that passes through the canopy without being absorbed or reflected. Deciduous trees in full leaf have transmittance values between 0.1 and 0.4 depending on species and LAI. 0.3 reflects moderate canopy density, consistent with empirical observations for mid-summer crowns. It controls how much light reaches the ground and sub-canopy vegetation; affects microclimate and shading.
Both values can be adjusted to match field measurements or literature for specific species or leaf phenology. However you can use them as robust fallback defaults when exact species traits are unavailable.
Season-Profile
Defines monthly LAD normalization.
SeasonProfile = c(0.2, 0.2, 0.4, 0.7, 1.0, 1.0, 1.0, 0.8, 0.6, 0.3, 0.2, 0.2)
The SeasonProfile is a vector of 12 numeric values (one per month) weighting the relative Leaf Area Density (LAD) throughout the year. It models seasonal leaf development and senescence, controlling how much foliage is present in each month:
- Values range from
0.0(no foliage) to1.0(full foliage). - For deciduous trees like Fagus sylvatica or Quercus robur, foliage develops in spring (April–May), peaks in summer (June–August), and declines in autumn (September–October).
Profile Breakdown:
| Months | Value | Interpretation |
|---|---|---|
| Jan–Feb, Nov–Dec | 0.2 | Dormant / leafless |
| March | 0.4 | Budburst begins |
| April | 0.7 | Leaf expansion |
| May–July | 1.0 | Full canopy |
| August | 0.8 | Leaf maturity decline |
| September | 0.6 | Senescence onset |
| October | 0.3 | Strong senescence |
The SeasonProfile directly influences LAD in ENVI-met’s dynamic vegetation simulation — affecting transpiration, shading, and energy balance across the simulation year. Adjusting this vector allows tailoring of phenology to site-specific or species-specific data.
L-SystemBased trees in ENVI-met (Experimetal)
ENVI-met optionally allows procedural generation of tree architecture using Lindenmayer Systems (L-Systems) — a formal grammar originally used to simulate plant growth patterns. When L-SystemBased = 1, the geometry of the tree is not derived from a static LAD matrix alone, but supplemented or replaced by rule-based 3D branching structures which supplement or replace the matrix. This is independent of the LAD profile but may affect shading and visualisation in the Albero interface of ENVI-met.
L-SystemBased = 1
Axiom = "F(2)V\V\\V/////B"
IterationDepth = 3Explanation of Key Parameters
| Parameter | Meaning |
|---|---|
L-SystemBased |
If 1, enables L-system generation (uses rules to grow plant structure) |
Axiom |
Starting string (“seed”) for the L-system; defines base growth |
IterationDepth |
How many times to apply production rules; higher means more detail |
TermLString |
Optional: Final symbol to be drawn/rendered (e.g. “L”) |
ApplyTermLString |
If 1, interprets the TermLString; otherwise, renders entire string |
Default Settings
<L-SystemBased>1</L-SystemBased>
<Axiom>F(2)V\V\\V/////B</Axiom>
<IterationDepth>3</IterationDepth>
<TermLString>L</TermLString>
<ApplyTermLString>1</ApplyTermLString>F(2): Move forward with length 2 (main trunk)V\\V/////B: Branching pattern with rotations (backslashes and slashes encode rotation commands);Bmay denote a terminal leaf or budIterationDepth = 3: The production rules (if defined) will be applied 3 times to this axiom, generating a fractal-like tree structure.
Note: In ENVI-met, the actual grammar rules are hard-coded and not customizable in
.pld— only the axiom and iteration depth are user-defined. It is highly experimental and poorly documented
Use L-SystemBased = 1 if:
You want visual structure added to otherwise sparse or low-resolution LAD matrices
The tree lacks realistic shape (for Albero visualization)
Use
L-SystemBased = 0(default) if:- You already provide a dense voxel-based LAD (from TLS or similar)
- You want strict control over the 3D structure via LAD profile only
Import TLS-based .pld into ENVI-met via Albero Clipboard
Requirements
- ENVI-met 5.8+
- .pld file (e.g. oak_tls_envimet.pld)
- Albero editor (via Leonardo)
Steps
1. Open Albero
→ Leonardo → Database → Plant Database
2. Open Clipboard
→ Click Clipboard (top-right)
3. Import .pld
→ Clipboard → Import → Load file
4. Edit (optional)
→ Adjust LAD, albedo, transmittance, name, etc.
5. Send to Library
→ Click “Send to Library”
6. Use in ENVI-met
→ In Leonardo/Spaces assign plant to your 3D model
Notes
- .pld contains LAD(z) values (m²/m³)
- Use Advanced Settings to fine-tune visualization
- Custom plants stored in your personal Albero library
Key Benefits
Efficient and scalable: The method avoids destructive sampling by using TLS return counts as proxies for leaf density. This makes it suitable for large-scale or repeated surveys without the need for time-consuming ground calibration.
Captures structural patterns: Normalizing the LAD values retains the vertical and spatial structure of vegetation, enabling meaningful comparison of crown shape, canopy layering, and vegetation density across space or time.
Directly usable in ENVI-met: The output is structured as a raster stack with height-specific layers, aligning with the input requirements of ENVI-met’s SimplePlant or 3D vegetation modules. This enables seamless integration into microclimate simulations.
Limitations
Simplified assumptions: The linear mapping of TLS returns to LAD assumes a proportional relationship, which simplifies the complex interaction between laser pulses and vegetation surfaces.
Scan geometry dependency: Occlusion, scan angle, and varying point densities can distort the return distribution, especially in dense or multi-layered vegetation.
Generic LAD normalization: The maximum LAD value used for normalization is taken from literature-based estimates rather than site-specific measurements, which can introduce bias in absolute LAD magnitudes.
Conclusion
This workflow offers a robust and accessible approach for analyzing vegetation structure and generating model-ready LAD profiles from TLS data. It is especially useful for relative comparisons and ecological modeling, but is not intended for absolute LAD quantification without additional calibration.
References
- Calders et al. (2015). Nondestructive biomass estimation via TLS. Methods Ecol Evol, 6:198–208.https://doi.org/10.1111/2041-210X.12301
- Chen et al. (2018): Estimation of LAI in open-canopy forests using TLS and path length models. Agric. For. Meteorol. 263, 323–333. https://doi.org/10.1016/j.agrformet.2018.09.006
- ENVI-met PLANT3D specification: https://www.envi-met.net/documents/papers/overview30.pdf
- ENVI-met Albero overview: https://envi-met.com/tutorials/albero-overview
- ENVI-met KB – Obtaining Leaf Area Density: https://envi-met.info/doku.php?id=kb:lad#obtaining_leaf_area_density_data
- ENVI-met dbmanager documentation: https://envi-met.info/doku.php?id=apps:dbmanager:start
- ENVI-met Vegetation Tutorial (YouTube): https://www.youtube.com/watch?v=KGRLnXAXZds
- Flynn et al. (2023) – TLS-based vegetation index estimation; compares methods and highlights complexities in Mediterranean forest. Biogeosciences, 20(13), 2769–2784. doi:10.5194/bg-20-2769-2023
- Hosoi & Omasa (2006). Voxel-based 3D tree modeling. IEEE TGRS, 44(12), 3610–3618. https://doi.org/10.1109/TGRS.2006.881743
- Prusinkiewicz & Lindenmayer (1990). The Algorithmic Beauty of Plants. Springer. https://doi.org/10.1007/978-1-4613-8476-2
- Oshio & Asawa (2016). Solar transmittance of urban trees. IEEE TGRS, 54(9), 5483–5492. https://doi.org/10.1109/TGRS.2016.2565699
- Simon, Sinsel & Bruse (2020). Fractal trees in ENVI-met. Forests, 11(8), 869. https://doi.org/10.3390/f11080869
- Wilkes et al. (2017). TLS acquisition strategies. Remote Sens Environ, 196, 140–153. https://doi.org/10.1016/j.rse.2017.04.030
- Chen et al. (2018). LAI from TLS. Agr Forest Meteorol, 263, 323–333. https://doi.org/10.1016/j.agrformet.2018.09.006
- Yin et al. (2019). Shading and thermal comfort. Sustainability, 11(5), 1355. https://doi.org/10.3390/su11051355
- Zhang (2024). Green layouts in ENVI-met. Informatica, 48(23). https://doi.org/10.31449/inf.v48i23.6881 Certainly. Here’s the reference adapted to match your current compact style: